{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "</center>\n",
    "Авторы материала: Data Science интерн Ciklum, студент магистерской программы CSDS UCU Виталий Радченко, программист-исследователь Mail.ru Group, старший преподаватель Факультета Компьютерных Наук ВШЭ Юрий Кашницкий. Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Тема 5. Композиции алгоритмов, случайный лес</center>\n",
    "## <center> Часть 2. Случайный лес</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лео Брейман нашел применение бутстрэпу не только в статистике, но и в машинном обучении. Он вместе с Адель Катлер усовершенстовал алгоритм случайного леса, который был предложенный [Хо](http://ect.bell-labs.com/who/tkh/publications/papers/odt.pdf), добавив к первоначальному варианту построение некоррелируемых деревьев на основе CART, в сочетании с методом случайных подпространств и бэггинга. \n",
    "\n",
    "Решающие деревья являются хорошим семейством базовых классификаторов для бэггинга, поскольку они достаточно сложны и могут достигать нулевой ошибки на любой выборке. Метод случайных подпространств позволяет снизить коррелированность между деревьями и избежать переобучения. Базовые алгоритмы обучаются на различных подмножествах признакового описания, которые также выделяются случайным образом.\n",
    "Ансамбль моделей, использующих метод случайного подпространства, можно построить, используя следующий алгоритм:\n",
    "1. Пусть количество объектов для тренировки равно $\\large N$, а количество признаков $\\large D$.\n",
    "2. Выберите $\\large L$ как число отдельных моделей в ансамбле.\n",
    "3. Для каждой отдельной модели $\\large l$ выберите $\\large dl (dl < D) $ как число признаков для $\\large l$ . Обычно для всех моделей используется только одно значение $\\large dl$.\n",
    "4. Для каждой отдельной модели $\\large l$ создайте обучающую выборку, выбрав $\\large dl$-признаков из $\\large D$ с заменой и обучите модель.\n",
    "\n",
    "Теперь, чтобы применить модель ансамбля к новому объекту, объедините результаты отдельных $\\large L$ моделей мажоритарным голосованием или путем комбинирования апостериорных вероятностей.\n",
    "\n",
    "\n",
    "### Алгоритм\n",
    "\n",
    "Алгоритм построения случайного леса, состоящего из $\\large N$ деревьев, выглядит следующим образом:\n",
    "* Для каждого $\\large n = 1, \\dots, N$:\n",
    "     * Сгенерировать выборку $\\large X_n$ с помощью bootstrap.\n",
    "     * Построить решающее дерево $\\large b_n$ по выборке $\\large X_n$:\n",
    "         — по заданному критерию мы выбираем лучший признак, делаем разбиение в дереве по нему и так до исчерпания выборки\n",
    "         — дерево строится, пока в каждом листе не более $\\large n_\\text{min}$ объектов или пока не достигнем определенной высоты дерева\n",
    "         — при каждом разбиении сначала выбирается $\\large m$ случайных признаков из $\\large n$ исходных, \n",
    "         и оптимальное разделение выборки ищется только среди них.\n",
    "         \n",
    "Итоговый классификатор $\\large a(x) = \\frac{1}{N}\\sum_{i = 1}^N b_i(x)$, простыми словами — для задачи кассификации мы выбираем решение голосованием по большинству, а в задаче регрессии — средним.\n",
    "\n",
    "Рекомендуется в задачах классификации брать $\\large m = \\sqrt{n}$, а в задачах регрессии — $\\large m = \\frac{n}{3}$, где $\\large n$ — число признаков. Также рекомендуется в задачах классификации строить каждое дерево до тех пор, пока в каждом листе не окажется по одному объекту, а в задачах регрессии — пока в каждом листе не окажется по пять объектов.\n",
    "\n",
    "Таким образом, случайный лес — это бэггинг над решающими деревьями, при обучении которых для каждого разбиения признаки выбираются из некоторого случайного подмножества признаков.\n",
    "\n",
    "### Сравнение с деревом решений и бэггингом"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "\n",
    "# отключим всякие предупреждения Anaconda\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = 10, 6\n",
    "import seaborn as sns\n",
    "from sklearn.datasets import make_circles\n",
    "from sklearn.ensemble import (BaggingClassifier, BaggingRegressor,\n",
    "                              RandomForestClassifier, RandomForestRegressor)\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor\n",
    "\n",
    "n_train = 150        \n",
    "n_test = 1000       \n",
    "noise = 0.1\n",
    "\n",
    "# Generate data\n",
    "def f(x):\n",
    "    x = x.ravel()\n",
    "    return np.exp(-x ** 2) + 1.5 * np.exp(-(x - 2) ** 2)\n",
    "\n",
    "def generate(n_samples, noise):\n",
    "    X = np.random.rand(n_samples) * 10 - 5\n",
    "    X = np.sort(X).ravel()\n",
    "    y = np.exp(-X ** 2) + 1.5 * np.exp(-(X - 2) ** 2)\\\n",
    "        + np.random.normal(0.0, noise, n_samples)\n",
    "    X = X.reshape((n_samples, 1))\n",
    "\n",
    "    return X, y\n",
    "\n",
    "X_train, y_train = generate(n_samples=n_train, noise=noise)\n",
    "X_test, y_test = generate(n_samples=n_test, noise=noise)\n",
    "\n",
    "# One decision tree regressor\n",
    "dtree = DecisionTreeRegressor().fit(X_train, y_train)\n",
    "d_predict = dtree.predict(X_test)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(X_test, f(X_test), \"b\")\n",
    "plt.scatter(X_train, y_train, c=\"b\", s=20)\n",
    "plt.plot(X_test, d_predict, \"g\", lw=2)\n",
    "plt.xlim([-5, 5])\n",
    "plt.title(\"Решающее дерево, MSE = %.2f\" \n",
    "          % np.sum((y_test - d_predict) ** 2))\n",
    "\n",
    "# Bagging decision tree regressor\n",
    "bdt = BaggingRegressor(DecisionTreeRegressor()).fit(X_train, y_train)\n",
    "bdt_predict = bdt.predict(X_test)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(X_test, f(X_test), \"b\")\n",
    "plt.scatter(X_train, y_train, c=\"b\", s=20)\n",
    "plt.plot(X_test, bdt_predict, \"y\", lw=2)\n",
    "plt.xlim([-5, 5])\n",
    "plt.title(\"Бэггинг решающих деревьев, MSE = %.2f\" % np.sum((y_test - bdt_predict) ** 2));\n",
    "\n",
    "# Random Forest\n",
    "rf = RandomForestRegressor(n_estimators=10).fit(X_train, y_train)\n",
    "rf_predict = rf.predict(X_test)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(X_test, f(X_test), \"b\")\n",
    "plt.scatter(X_train, y_train, c=\"b\", s=20)\n",
    "plt.plot(X_test, rf_predict, \"r\", lw=2)\n",
    "plt.xlim([-5, 5])\n",
    "plt.title(\"Случайный лес, MSE = %.2f\" % np.sum((y_test - rf_predict) ** 2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как мы видим из графиков и значений ошибки MSE, случайный лес из 10 деревьев дает лучший результат, чем одно дерево или бэггинг из 10 деревьев решений. Основное различие случайного леса и бэггинга на деревьях решений заключается в том, что в случайном лесе из случайного числа выбирается случайное подмножество признаков, и лучший признак для разделения узла определяется из подвыборки признаков, в отличие от бэггинга, где все функции рассматриваются для разделения в узле.\n",
    "\n",
    "Также можно увидеть преимущество случайного леса и бэггинга в задачах классификации.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "X, y = make_circles(n_samples=500, factor=0.1, noise=0.35, random_state=42)\n",
    "X_train_circles, X_test_circles, y_train_circles, y_test_circles = train_test_split(X, y, test_size=0.2)\n",
    "\n",
    "dtree = DecisionTreeClassifier(random_state=42)\n",
    "dtree.fit(X_train_circles, y_train_circles)\n",
    "\n",
    "x_range = np.linspace(X.min(), X.max(), 100)\n",
    "xx1, xx2 = np.meshgrid(x_range, x_range)\n",
    "y_hat = dtree.predict(np.c_[xx1.ravel(), xx2.ravel()])\n",
    "y_hat = y_hat.reshape(xx1.shape)\n",
    "plt.contourf(xx1, xx2, y_hat, alpha=0.2)\n",
    "plt.scatter(X[:,0], X[:,1], c=y, cmap='autumn')\n",
    "plt.title(\"Дерево решений\")\n",
    "plt.show()\n",
    "\n",
    "b_dtree = BaggingClassifier(DecisionTreeClassifier(),n_estimators=300, random_state=42)\n",
    "b_dtree.fit(X_train_circles, y_train_circles)\n",
    "\n",
    "x_range = np.linspace(X.min(), X.max(), 100)\n",
    "xx1, xx2 = np.meshgrid(x_range, x_range)\n",
    "y_hat = b_dtree.predict(np.c_[xx1.ravel(), xx2.ravel()])\n",
    "y_hat = y_hat.reshape(xx1.shape)\n",
    "plt.contourf(xx1, xx2, y_hat, alpha=0.2)\n",
    "plt.scatter(X[:,0], X[:,1], c=y, cmap='autumn')\n",
    "plt.title(\"Бэггинг(дерево решений)\")\n",
    "plt.show()\n",
    "\n",
    "rf = RandomForestClassifier(n_estimators=300, random_state=42)\n",
    "rf.fit(X_train_circles, y_train_circles)\n",
    "\n",
    "x_range = np.linspace(X.min(), X.max(), 100)\n",
    "xx1, xx2 = np.meshgrid(x_range, x_range)\n",
    "y_hat = rf.predict(np.c_[xx1.ravel(), xx2.ravel()])\n",
    "y_hat = y_hat.reshape(xx1.shape)\n",
    "plt.contourf(xx1, xx2, y_hat, alpha=0.2)\n",
    "plt.scatter(X[:,0], X[:,1], c=y, cmap='autumn')\n",
    "plt.title(\"Случайный лес\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "На рисунках выше видно, что разделяющая граница дерева решений очень «рваная» и на ней много острых углов, что говорит о переобучении и слабой обобщающей способности. В то время как у бэггинга и случайного леса граница достаточно сглаженная и практически нет признаков переобучения.\n",
    "\n",
    "Давайте теперь попробуем разобраться с параметрами, с помощью подбора которых мы cможем увеличить долю правильных ответов.\n",
    "\n",
    "### Параметры\n",
    "\n",
    "\n",
    "Метод случайного леса реализован в библиотеке машинного обучения [scikit-learn](http://scikit-learn.org/stable/) двумя классами RandomForestClassifier и RandomForestRegressor.\n",
    "\n",
    "Полный список параметров случайного леса для задачи регрессии:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```class sklearn.ensemble.RandomForestRegressor(\n",
    "    n_estimators — число деревьев в \"лесу\" (по дефолту – 10)\n",
    "    criterion — функция, которая измеряет качество разбиения ветки дерева (по дефолту — \"mse\" , так же можно выбрать \"mae\")\n",
    "    max_features — число признаков, по которым ищется разбиение. Вы можете указать конкретное число или процент признаков, либо выбрать из доступных значений: \"auto\" (все признаки), \"sqrt\", \"log2\". По дефолту стоит \"auto\".\n",
    "    max_depth — максимальная глубина дерева  (по дефолту глубина не ограничена)\n",
    "    min_samples_split — минимальное количество объектов, необходимое для разделения внутреннего узла. Можно задать числом или процентом от общего числа объектов (по дефолту — 2)\n",
    "    min_samples_leaf — минимальное число объектов в листе. Можно задать числом или процентом от общего числа объектов (по дефолту — 1)\n",
    "    min_weight_fraction_leaf — минимальная взвешенная доля от общей суммы весов (всех входных объектов) должна быть в листе (по дефолту имеют одинаковый вес)\n",
    "    max_leaf_nodes — максимальное количество листьев (по дефолту нет ограничения)\n",
    "    min_impurity_split — порог для остановки наращивания дерева (по дефолту 1е-7)\n",
    "    bootstrap — применять ли бустрэп для построения дерева (по дефолту True)\n",
    "    oob_score — использовать ли out-of-bag объекты для оценки R^2 (по дефолту False)\n",
    "    n_jobs — количество ядер для построения модели и предсказаний (по дефолту 1, если поставить -1, то будут использоваться все ядра)\n",
    "    random_state — начальное значение для генерации случайных чисел (по дефолту его нет, если хотите воспроизводимые результаты, то нужно указать любое число типа int\n",
    "    verbose — вывод логов по построению деревьев (по дефолту 0)\n",
    "    warm_start — использует уже натренированую модель и добавляет деревьев в ансамбль (по дефолту False)\n",
    ")```\n",
    "\n",
    "Для задачи классификации все почти то же самое, мы приведем только те параметры, которыми RandomForestClassifier отличается от RandomForestRegressor \n",
    "\n",
    "```class sklearn.ensemble.RandomForestClassifier(\n",
    "    criterion — поскольку у нас теперь задача классификации, то по дефолту выбран критерий \"gini\" (можно выбрать \"entropy\")\n",
    "    class_weight — вес каждого класса (по дефолту все веса равны 1, но можно передать словарь с весами, либо явно указать \"balanced\", тогда веса классов будут равны их исходным частям в генеральной совокупности; также можно указать \"balanced_subsample\", тогда веса на каждой подвыборке будут меняться в зависимости от распределения классов на этой подвыборке.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Далее рассмотрим несколько параметров, на которые в первую очередь стоит обратить внимание при построении модели:\n",
    "- n_estimators — число деревьев в \"лесу\"\n",
    "- criterion — критерий для разбиения выборки в вершине\n",
    "- max_features — число признаков, по которым ищется разбиение\n",
    "- min_samples_leaf — минимальное число объектов в листе\n",
    "- max_depth — максимальная глубина дерева\n",
    "\n",
    "**Рассмотрим применение случайного леса в реальной задаче**\n",
    "\n",
    "Для этого будем использовать пример с задачей fraud detection. Это задача классификации, поэтому будем использовать метрику accuracy для оценки точности. Для начала построим самый простой классификатор, который будет нашим бейслайном. Возьмем только числовые признаки для упрощения.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.model_selection import (GridSearchCV, StratifiedKFold,\n",
    "                                     cross_val_score)\n",
    "\n",
    "# Загружаем данные\n",
    "df = pd.read_csv(\"../../data/telecom_churn.csv\")\n",
    "\n",
    "# Выбираем сначала только колонки с числовым типом данных\n",
    "cols = []\n",
    "for i in df.columns:\n",
    "    if (df[i].dtype == \"float64\") or (df[i].dtype == 'int64'):\n",
    "        cols.append(i)\n",
    "        \n",
    "# Разделяем на признаки и объекты\n",
    "X, y = df[cols].copy(), np.asarray(df[\"Churn\"],dtype='int8')\n",
    "\n",
    "# Инициализируем стратифицированную разбивку нашего датасета для валидации\n",
    "skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "# Инициализируем наш классификатор с дефолтными параметрами\n",
    "rfc = RandomForestClassifier(random_state=42, n_jobs=-1, oob_score=True)\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "results = cross_val_score(rfc, X, y, cv=skf)\n",
    "\n",
    "# Оцениваем точность на тестовом датасете\n",
    "print(\"CV accuracy score: {:.2f}%\".format(results.mean()*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Получили точность 91.21%, теперь попробуем улучшить этот результат и посмотреть, как ведут себя кривые обучения при изменении основных параметров.\n",
    "\n",
    "Начнем с количества деревьев:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Инициализируем валидацию\n",
    "skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n",
    "\n",
    "# Создаем списки для сохранения точности на тренировочном и тестовом датасете\n",
    "train_acc = []\n",
    "test_acc = []\n",
    "temp_train_acc = []\n",
    "temp_test_acc = []\n",
    "trees_grid = [5, 10, 15, 20, 30, 50, 75, 100]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for ntrees in trees_grid:\n",
    "    rfc = RandomForestClassifier(n_estimators=ntrees, random_state=42, n_jobs=-1, oob_score=True)\n",
    "    temp_train_acc = []\n",
    "    temp_test_acc = []\n",
    "    for train_index, test_index in skf.split(X, y):\n",
    "        X_train, X_test = X.iloc[train_index], X.iloc[test_index]\n",
    "        y_train, y_test = y[train_index], y[test_index]\n",
    "        rfc.fit(X_train, y_train)\n",
    "        temp_train_acc.append(rfc.score(X_train, y_train))\n",
    "        temp_test_acc.append(rfc.score(X_test, y_test))\n",
    "    train_acc.append(temp_train_acc)\n",
    "    test_acc.append(temp_test_acc)\n",
    "    \n",
    "train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)\n",
    "print(\"Best accuracy on CV is {:.2f}% with {} trees\".format(max(test_acc.mean(axis=1))*100, \n",
    "                                                        trees_grid[np.argmax(test_acc.mean(axis=1))]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.style.use('ggplot')\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(trees_grid, train_acc.mean(axis=1), alpha=0.5, color='blue', label='train')\n",
    "ax.plot(trees_grid, test_acc.mean(axis=1), alpha=0.5, color='red', label='cv')\n",
    "ax.fill_between(trees_grid, test_acc.mean(axis=1) - test_acc.std(axis=1), test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)\n",
    "ax.fill_between(trees_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1), test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)\n",
    "ax.legend(loc='best')\n",
    "ax.set_ylim([0.88,1.02])\n",
    "ax.set_ylabel(\"Accuracy\")\n",
    "ax.set_xlabel(\"N_estimators\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видно, при достижении определенного числа деревьев, наша точность на тесте выходит на асимптоту, и вы можете сами решить, сколько деревьев оптимально для вашей задачи. \n",
    "На рисунке также видно, что на тренировочной выборке мы смогли достичь 100% точности, это говорит нам о переобучении нашей модели. Чтобы избежать переобучения, мы должны добавить параметры регуляризации в модель.\n",
    "\n",
    "Начнем с параметра максимальной глубины – `max_depth`. (зафиксируем к-во деревьев 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Создаем списки для сохранения точности на тренировочном и тестовом датасете\n",
    "train_acc = []\n",
    "test_acc = []\n",
    "temp_train_acc = []\n",
    "temp_test_acc = []\n",
    "max_depth_grid = [3, 5, 7, 9, 11, 13, 15, 17, 20, 22, 24]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for max_depth in max_depth_grid:\n",
    "    rfc = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1, oob_score=True, max_depth=max_depth)\n",
    "    temp_train_acc = []\n",
    "    temp_test_acc = []\n",
    "    for train_index, test_index in skf.split(X, y):\n",
    "        X_train, X_test = X.iloc[train_index], X.iloc[test_index]\n",
    "        y_train, y_test = y[train_index], y[test_index]\n",
    "        rfc.fit(X_train, y_train)\n",
    "        temp_train_acc.append(rfc.score(X_train, y_train))\n",
    "        temp_test_acc.append(rfc.score(X_test, y_test))\n",
    "    train_acc.append(temp_train_acc)\n",
    "    test_acc.append(temp_test_acc)\n",
    "    \n",
    "train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)\n",
    "print(\"Best accuracy on CV is {:.2f}% with {} max_depth\".format(max(test_acc.mean(axis=1))*100, \n",
    "                                                        max_depth_grid[np.argmax(test_acc.mean(axis=1))]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(max_depth_grid, train_acc.mean(axis=1), alpha=0.5, color='blue', label='train')\n",
    "ax.plot(max_depth_grid, test_acc.mean(axis=1), alpha=0.5, color='red', label='cv')\n",
    "ax.fill_between(max_depth_grid, test_acc.mean(axis=1) - test_acc.std(axis=1), test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)\n",
    "ax.fill_between(max_depth_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1), test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)\n",
    "ax.legend(loc='best')\n",
    "ax.set_ylim([0.88,1.02])\n",
    "ax.set_ylabel(\"Accuracy\")\n",
    "ax.set_xlabel(\"Max_depth\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Параметр `max_depth` хорошо справляется с регуляризацией модели, и мы уже не так сильно переобучаемся. Точность нашей модели немного возросла.\n",
    "\n",
    "Еще важный параметр `min_samples_leaf`, он так же выполняет функцию регуляризатора.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Создаем списки для сохранения точности на тренировочном и тестовом датасете\n",
    "train_acc = []\n",
    "test_acc = []\n",
    "temp_train_acc = []\n",
    "temp_test_acc = []\n",
    "min_samples_leaf_grid = [1, 3, 5, 7, 9, 11, 13, 15, 17, 20, 22, 24]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for min_samples_leaf in min_samples_leaf_grid:\n",
    "    rfc = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1, \n",
    "                                 oob_score=True, min_samples_leaf=min_samples_leaf)\n",
    "    temp_train_acc = []\n",
    "    temp_test_acc = []\n",
    "    for train_index, test_index in skf.split(X, y):\n",
    "        X_train, X_test = X.iloc[train_index], X.iloc[test_index]\n",
    "        y_train, y_test = y[train_index], y[test_index]\n",
    "        rfc.fit(X_train, y_train)\n",
    "        temp_train_acc.append(rfc.score(X_train, y_train))\n",
    "        temp_test_acc.append(rfc.score(X_test, y_test))\n",
    "    train_acc.append(temp_train_acc)\n",
    "    test_acc.append(temp_test_acc)\n",
    "    \n",
    "train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)\n",
    "print(\"Best accuracy on CV is {:.2f}% with {} min_samples_leaf\".format(max(test_acc.mean(axis=1))*100, \n",
    "                                                        min_samples_leaf_grid[np.argmax(test_acc.mean(axis=1))]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(min_samples_leaf_grid, train_acc.mean(axis=1), alpha=0.5, color='blue', label='train')\n",
    "ax.plot(min_samples_leaf_grid, test_acc.mean(axis=1), alpha=0.5, color='red', label='cv')\n",
    "ax.fill_between(min_samples_leaf_grid, test_acc.mean(axis=1) - test_acc.std(axis=1), test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)\n",
    "ax.fill_between(min_samples_leaf_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1), test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)\n",
    "ax.legend(loc='best')\n",
    "ax.set_ylim([0.88,1.02])\n",
    "ax.set_ylabel(\"Accuracy\")\n",
    "ax.set_xlabel(\"Min_samples_leaf\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В данном случае мы не выигрываем в точности на валидации, но зато можем сильно уменьшить переобучение до 2% при сохранении точности около 92%.\n",
    "\n",
    "Рассмотрим такой параметр как `max_features`. Для задач классификации по умолчанию используется $\\large \\sqrt{n}$, где n — число признаков. Давайте проверим, оптимально ли в нашем случае использовать 4 признака или нет.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Создаем списки для сохранения точности на тренировочном и тестовом датасете\n",
    "train_acc = []\n",
    "test_acc = []\n",
    "temp_train_acc = []\n",
    "temp_test_acc = []\n",
    "max_features_grid = [2, 4, 6, 8, 10, 12, 14, 16]\n",
    "\n",
    "# Обучаем на тренировочном датасете\n",
    "for max_features in max_features_grid:\n",
    "    rfc = RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1, \n",
    "                                 oob_score=True, max_features=max_features)\n",
    "    temp_train_acc = []\n",
    "    temp_test_acc = []\n",
    "    for train_index, test_index in skf.split(X, y):\n",
    "        X_train, X_test = X.iloc[train_index], X.iloc[test_index]\n",
    "        y_train, y_test = y[train_index], y[test_index]\n",
    "        rfc.fit(X_train, y_train)\n",
    "        temp_train_acc.append(rfc.score(X_train, y_train))\n",
    "        temp_test_acc.append(rfc.score(X_test, y_test))\n",
    "    train_acc.append(temp_train_acc)\n",
    "    test_acc.append(temp_test_acc)\n",
    "    \n",
    "train_acc, test_acc = np.asarray(train_acc), np.asarray(test_acc)\n",
    "print(\"Best accuracy on CV is {:.2f}% with {} max_features\".format(max(test_acc.mean(axis=1))*100, \n",
    "                                                        max_features_grid[np.argmax(test_acc.mean(axis=1))]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "ax.plot(max_features_grid, train_acc.mean(axis=1), alpha=0.5, color='blue', label='train')\n",
    "ax.plot(max_features_grid, test_acc.mean(axis=1), alpha=0.5, color='red', label='cv')\n",
    "ax.fill_between(max_features_grid, test_acc.mean(axis=1) - test_acc.std(axis=1), test_acc.mean(axis=1) + test_acc.std(axis=1), color='#888888', alpha=0.4)\n",
    "ax.fill_between(max_features_grid, test_acc.mean(axis=1) - 2*test_acc.std(axis=1), test_acc.mean(axis=1) + 2*test_acc.std(axis=1), color='#888888', alpha=0.2)\n",
    "ax.legend(loc='best')\n",
    "ax.set_ylim([0.88,1.02])\n",
    "ax.set_ylabel(\"Accuracy\")\n",
    "ax.set_xlabel(\"Max_features\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В нашем случае оптимальное число признаков — 10, именно с таким значением достигается наилучший результат.\n",
    "\n",
    "Мы рассмотрели, как ведут себя кривые обучения в зависимости от изменения основных параметров. Давайте теперь с помощью `GridSearch` найдем оптимальные параметры для нашего примера.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Сделаем инициализацию параметров, по которым хотим сделать полный перебор\n",
    "parameters = {'max_features': [4, 7, 10, 13], 'min_samples_leaf': [1, 3, 5, 7], 'max_depth': [5,10,15,20]}\n",
    "rfc = RandomForestClassifier(n_estimators=100, random_state=42, \n",
    "                             n_jobs=-1, oob_score=True)\n",
    "gcv = GridSearchCV(rfc, parameters, n_jobs=-1, cv=skf, verbose=1)\n",
    "gcv.fit(X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcv.best_estimator_, gcv.best_score_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Лучшая точность, который мы смогли достичь с помощью перебора параметров — 92.83% при `'max_depth': 15, 'max_features': 7, 'min_samples_leaf': 3`. \n",
    "\n",
    "### Вариация и декорреляционный эффект\n",
    "\n",
    "Давайте запишем вариацию для случайного леса как $$ \\large Varf(x) = \\rho(x)\\sigma^2(x) $$\n",
    "Тут  \n",
    "- $ \\large \\rho(x)$ – корреляция выборки между любыми двумя деревьями, используемыми при усреднении \n",
    "$$ \\large \\rho(x) = corr[T(x;\\Theta_1(Z)),T(x_2,\\Theta_2(Z))], $$ где $ \\large \\Theta_1(Z) $ и $ \\large  \\Theta_2(Z) $ – случайно выбранная пара деревьев на случайно выбранных объектах выборки Z\n",
    "- $ \\large \\sigma^2(x)$ — это выборочная дисперсия любого произвольно выбранного дерева: $$ \\large \\sigma^2(x) = VarT(x;\\Theta(X))$$\n",
    "\n",
    "Легко спутать $ \\large \\rho(x) $ со средней корреляцией между обученными деревьями в данном случайном лесе, рассматривая деревья как N-векторы и вычисляя среднюю парную корреляцию между ними. Это не тот случай. Эта условная корреляция не имеет прямого отношения к процессу усреднения, а зависимость от  $ \\large x$  в $ \\large \\rho(x)$ предупреждает нас об этом различии. Скорее, $ \\large \\rho(x)$ является теоретической корреляцией между парой случайных деревьев, оцененных в объекте $ \\large x$, которая была вызвана многократным сэмплированием тренировочной выборки из генеральной совокупности $ \\large Z $, и после этого выбрана данная пара случайных деревьев. На статистическом жаргоне это корреляция, вызванная выборочным распределением $ \\large Z $ и $ \\large \\Theta $.\n",
    "\n",
    "По факту, условная ковариация пары деревьев равна 0, потому что бустрэп и отбор признаков — независимы и одинаково распределены.\n",
    "\n",
    "Если рассмотреть дисперсию по одному дереву, то она практически не меняется от переменных для разделения ($ \\large m $), а вот для ансамбля это играет большую роль, и дисперсия для дерева намного выше, чем для ансамбля. \n",
    "В книге *The Elements of Statistical Learning (Trevor Hastie, Robert Tibshirani и Jerome Friedman)* есть отличный пример, который это демонстрирует.\n",
    "![image](../../img/variance_rf.png)\n",
    "\n",
    "### Смещение\n",
    "\n",
    "Как и в бэггинге, смещение в случайном лесе такое же, как и смещение в отдельно взятом дереве $ \\large T(x,\\Theta(Z))$ \n",
    "$$ \\large \\begin{array}{rcl} Bias &=& \\mu(x) - E_Zf_{rf}(x) \\\\\n",
    "&=& \\mu(x) - E_ZE_{\\Theta | Z}T(x,\\Theta(Z))\\end{array}$$\n",
    "Это также обычно больше (в абсолютных величинах), чем смещение «неусеченного» (unprunned) дерева, поскольку рандомизация и сокращение пространства выборки налагают ограничения. Следовательно, улучшения в прогнозировании, полученные с помощью бэггинга или случайных лесов, являются исключительно результатом уменьшения дисперсии.\n",
    "\n",
    "###  Сверхслучайные деревья\n",
    "\n",
    "В сверхслучайных деревьях (Extremely Randomized Trees) больше случайности в том, как вычисляются разделения в узлах. Как и в случайных лесах, используется случайное подмножество возможных признаков, но вместо поиска наиболее оптимальных порогов, пороговые значения произвольно выбираются для каждого возможного признака, и наилучший из этих случайно генерируемых порогов выбирается как лучшее правило для разделения узла. Это обычно позволяет немного уменьшить дисперсию модели за счет несколько большего увеличения смещения.\n",
    "\n",
    "В библиотеке scikit-learn есть реализация [ExtraTreesClassifier](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.ExtraTreesClassifier.html#sklearn.ensemble.ExtraTreesClassifier) и [ExtraTreesRegressor](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.ExtraTreesRegressor.html#sklearn.ensemble.ExtraTreesRegressor). Данный метод стоит использовать, когда вы сильно переобучаетесь на случайном лесе или градиентном бустинге. \n",
    "\n",
    "### Схожесть случайного леса с алгоритмом k-ближайших соседей**\n",
    "\n",
    "Метод случайного леса схож с методом ближайших соседей. Случайные леса, по сути, осуществляют предсказания для объектов на основе меток похожих объектов из обучения. Схожесть объектов при этом тем выше, чем чаще эти объекты оказываются в одном и том же листе дерева. Покажем это формально.\n",
    "\n",
    "Рассмотрим задачу регрессии с квадратичной функцией потерь. Пусть $ \\large T_n(x) $ — номер листа n-го дерева из случайного леса, в который попадает объект $ \\large x $ . Ответ объекта $ \\large x $ равен среднему ответу по всем объектам обучающей выборки, которые попали в этот лист $ \\large T_n(x) $. Это можно записать так \n",
    "$$ \\large b_n(x) = \\sum_{i=1}^{l}w_n(x,x_i)y_i,$$ \n",
    "где $$ \\large w_n(x, x_i) = \\frac{[T_n(x) = T_n(x_i)]}{\\sum_{j=1}^{l}[T_n(x) = T_n(x_j)]}$$\n",
    "Тогда ответ композиции равен $$ \\large \\begin{array}{rcl} a_n(x) &=& \\frac{1}{N}\\sum_{n=1}^{N}\\sum_{i=1}^{l}w_n(x,x_i)y_i \\\\\n",
    "&=& \\sum_{i=1}^{l}\\Big(\\frac{1}{N}\\sum_{j=1}^{N}w_n(x,x_j)\\Big)y_i \\end{array}$$\n",
    "Видно, что ответ случайного леса представляет собой сумму ответов всех объектов обучения с некоторыми весами. Отметим, что номер листа $ \\large T_n(x)$, в который попал объект, сам по себе является ценным признаком. Достаточно неплохо работает подход, в котором по выборке обучается композиция из небольшого числа деревьев с помощью случайного леса или градиентного бустинга, а потом к ней добавляются категориальные признаки $ \\large T_1(x), \\dots, T_n(x) $. Новые признаки являются результатом нелинейного разбиения пространства и несут в себе информацию о сходстве объектов.\n",
    "\n",
    "Все в той же книге *The Elements of Statistical Learning* есть хороший наглядный пример сходства случайного леса и k-ближайших соседей. \n",
    "\n",
    "![image](../../img/knn_vs_rf.png)\n",
    "\n",
    "### Преобразование признаков в многомерное пространство\n",
    "\n",
    "Все привыкли использовать случайный лес для задач обучения с учителем, но также есть возможность проводить обучение и без учителя. С помощью метода [RandomTreesEmbedding](http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomTreesEmbedding.html#sklearn.ensemble.RandomTreesEmbedding) мы можем сделать трансформацию нашего датасета в многомерное разреженное его представление. Его суть в том, что мы строим абсолютно случайные деревья, и индекс листа, в котором оказалось наблюдение, мы считаем за новый признак. Если в первый лист попал объект, то мы ставим 1, а если не попал, то 0. Так называемое бинарное кодирование. Контролировать количество переменных и также степень разряженности нашего нового представления датасета мы можем увеличивая/уменьшая количество деревьев и их глубины. Поскольку соседние точки данных, скорее всего, лежат в одном и том же листе дерева, преобразование выполняет неявную, непараметрическую оценку плотности.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Плюсы**:\n",
    " - имеет высокую точность предсказания, на большинстве задач будет лучше линейных алгоритмов; точность сравнима с точностью бустинга\n",
    " - практически не чувствителен к выбросам в данных из-за случайного сэмлирования\n",
    " - нечувствителен к масштабированию (и вообще к любым монотонным преобразованиям) значений признаков, связано с выбором случайных подпространств\n",
    " - не требует тщательной настройки параметров, хорошо работает «из коробки». С помощью «тюнинга» параметров можно достичь прироста от 0.5 до 3% точности в зависимости от задачи и данных\n",
    " - способен эффективно обрабатывать данные с большим числом признаков и классов\n",
    " - одинаково хорошо обрабатывет как непрерывные, так и дискретные признаки\n",
    " - редко переобучается, на практике добавление деревьев почти всегда только улучшает композицию, но на валидации, после достижения определенного количества деревьев, кривая обучения выходит на асимптоту\n",
    " - для случайного леса существуют методы оценивания значимости отдельных признаков в модели\n",
    " - хорошо работает с пропущенными данными; сохраняет хорошую точность, если большая часть данных пропущенна\n",
    " - предполагает возможность сбалансировать вес каждого класса на всей выборке, либо на подвыборке каждого дерева\n",
    " - вычисляет близость между парами объектов, которые могут использоваться при кластеризации, обнаружении выбросов или (путем масштабирования) дают интересные представления данных\n",
    " - возможности, описанные выше, могут быть расширены до немаркированных данных, что приводит к возможноти делать кластеризацию и визуализацию данных, обнаруживать выбросы\n",
    " - высокая параллелизуемость и масштабируемость.\n",
    " \n",
    "**Минусы**:\n",
    " - в отличие от одного дерева, результаты случайного леса сложнее интерпретировать\n",
    " - нет формальных выводов (p-values), доступных для оценки важности переменных\n",
    " - алгоритм работает хуже многих линейных методов, когда в выборке очень много разреженных признаков (тексты, Bag of words)\n",
    " - случайный лес не умеет экстраполировать, в отличие от той же линейной регрессии (но это можно считать и плюсом, так как не будет экстремальных значений в случае попадания выброса)\n",
    " - алгоритм склонен к переобучению на некоторых задачах, особенно на зашумленных данных\n",
    " - для данных, включающих категориальные переменные с различным количеством уровней, случайные леса предвзяты в пользу признаков с большим количеством уровней: когда у признака много уровней, дерево будет сильнее подстраиваться именно под эти признаки, так как на них можно получить более высокую точность\n",
    " - если данные содержат группы коррелированных признаков, имеющих схожую значимость для меток, то предпочтение отдается небольшим группам перед большими\n",
    " - больший размер получающихся моделей. Требуется $O(NK)$ памяти для хранения модели, где $K$ — число деревьев.\n",
    "\n",
    "**Полезные источники:**\n",
    "- 15 раздел книги “[Elements of Statistical Learning](https://web.stanford.edu/~hastie/Papers/ESLII.pdf)” Jerome H. Friedman, Robert Tibshirani, and Trevor Hastie\n",
    "- [Блог](https://alexanderdyakonov.wordpress.com/2016/11/14/случайный-лес-random-forest/) Александра Дьяконова\n",
    "- больше про практические применение случайного леса и других алгоритмов композиций в официальной документации [scikit-learn](http://scikit-learn.org/stable/modules/ensemble.html)\n",
    "- [Курс](https://github.com/esokolov/ml-course-hse) Евгения Соколова по машинному обучению (материалы на GitHub). Есть дополнительные практические задания для углубления ваших знаний."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
